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It is known that there is no phase transition down to zero temperature in the antiferromagnetic 
Ising model on spatially anisotropic triangular lattices, in which the exchange coupling of one 
direction is stronger than those of other two directions. In the model, the low-temperature physics 
is governed by domain-wall excitations (defects) residing on bonds of the strong-coupling direction. 
In this paper, we show that an additional small attractive interaction between defects leads to a 
Berezinskii-Kosterlitz-Thouless (BKT) transition at a finite temperature, by applying the Monte 
Carlo simulation. The BKT phase can be viewed as the phase with a quasi long-range order of 
defects. We determine the phase diagram in a wide parameter regime and argue the phase structure 
from statistical-mechanics and field-theory viewpoints. 
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Introduction.— Antiferromagnets on triangular lat- 
tices [1] have long been studied as representative systems 
with geometric frustration. 0, Q In frustrated magnets, 
spins cannot take their lowest-energy configuration and 
hence macroscopic many-body states are almost degen- 
erated at sufficiently low temperatures. This is a main 
reason why various frustrated magnets often lead to in- 
triguing phenomena. For continuous or quantum spin 
systems on triangular lattices, spin chiralities [4] such 
as S-p x S r > have offered interesting topics. In addition 
to chiralities, recently multiferroic properties [H, spin 
nematic orders and spin-liquid states [3] have been 
vividly explored as novel many-body states or phenom- 
ena. 

Triangular antiferromagnets made from discrete spin 
degrees of freedom (Ising, clock, Potts models, etc.) [H 
are also interesting frustrated systems. They are also 
useful as simple models of, for instance, charge orderings, 
electric dipoles, and crystal growth. For these systems, 
their macroscopically degenerated ground states can be 
often counted out, but small thermal or quantum fluctua- 
tions in the ground-state manifold can induce unexpected 
phenomena. Hence effects of those fluctuations have of- 
fered an attractive research field. In this paper, we focus 
on antiferromagnetic (AF) spin-^ Ising systems on trian- 
gular lattices. We start with the following Hamiltonian 
with a spatial anisotropy: 



H 



J1-J2 



Ji 



S r S r f -\- J2 ^ ^ S r S r ' , (1) 



where S r is the Ising spin on site r = (r x ,r y ) and take 
values ±1, and J\ > and J 2 > respectively de- 
note the AF exchange coupling constants along one direc- 
tion of the triangular lattice and along other two direc- 
tions. This model is illustrated in Fig.[TJa). As shown in 
the figure, the symbols (r,r') x and {r,r') y respectively 
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FIG. 1: (a) AF Ising model 7ij 1 -j 2 Q on triangular lattice, 
in which white and black circles respectively denote up and 
down spins (S r = 1 and -1). Red thick bonds show the po- 
sitions of defects, (b) Three types of defects (bl)-(b3) in the 
model (pQ). Values Ji — 2J2, Ji, and J\ + 2J2 are the en- 
ergies of defects measured from the state without any bond, 
(c) Modified Ising model JLj 1 -j 2 -j 3 ([2]) with an additional 
next-nearest-neighbor coupling J3. 



stand for nearest-neighbor Ising spin pair along the x 
and y directions. The ratio J2/J1 measures the spatial 
anisotropy. The point J2/J1 = corresponds to decou- 
pled AF-Ji Ising chains, while J2/J1 = 1 is the spatially- 
isotropic triangular Ising model. We will consider the 
wide range < J2/J1 < 1 hereafter. 

The study of the spatially anisotropic Ising model (pQ) 
has a long history [8[. Wannier has studied thermody- 
namic properties of the isotropic model (pQ) with J x = J 2 , 
by using transfer matrix method [9j. He has accurately 
computed the residual entropy for the macroscopically 
degenerated ground states, and shown that the internal 
energy has no singular point as a function of temperature 
T. For J\ > J 2 > 0, the ground states consist of all of 
the states where Ising chains along the y direction are 
independently Neel ordered. In fact, when Neel ordered 
chains align to the y direction, the inter-chain energy is 
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canceled out. As a result, the ground states are semi- 
macroscopically degenerated [lOj. Houtappel has exam- 
ined the general anisotropic case with J\ ^ J2, and de- 
rived the exact expression of the partition function [Til ] . 
Recently, Hotta, et. al. 12, [l3| have studied thermo- 



dynamic properties of the model (pp) with J2/J1 < 1, by 
using the exact solution of Ref. EH and Monte Carlo (MC) 
simulation. They have confirmed that there is no phase 
transition in the plane of T and J2/ J\{< 1)- They have 
pointed out that the low-energy physics of the anisotropic 
model (pQ) with J2/J1 < 1 are governed by domain- wall 
excitations (defects) as shown in Fig. QJa) and (b). The 
defects are the ferromagnetic bonds cutting Neel ordered 
state of each Ising chain, and they are classified into 
three kinds in the energetic sense, as in Fig. QJb); Good, 
bad, and worst defects. In a sufficiently low-temperature 
regime, bad and worst defects disappear and only good 
defects survive. The inter-chain coupling J 2 develops a 
short-range correlation between defects along the x di- 
rection, but it does not induce any phase transition. 

These results make us infer that defects form a order- 
ing pattern (crystallization of defects) as in electrons of 
Wigner crystals if we introduce a small interaction be- 
tween defects in the anisotropic model (pp) with J2/J1 < 
1. Following this expectation, we will study the possibil- 
ity of defect crystallization in this paper. To this end, we 
introduce an additional ferromagnetic interaction J3 < 
to the J1-J2 model (pQ). The new Hamiltonian is given 

by 



J\ — Jl — J3 



n 



J1-J2 



(2) 



The J3 term is a next-nearest-neighbor interaction along 
the x direction and is illustrated in Fig. QJc). The panel 
(c) shows that a ferromagnetic J3 term plays a role of an 
attraction between two neighboring good defects. 

Numerical analysis.— We will apply the standard heat- 
bath MC method to the J1-J2-J3 model ([2]). Both lengths 
of the x and y directions, L x and L yi are set to L, and 
the periodic boundary condition for both directions are 
imposed. The system size L is increased up to 128. Tak- 
ing the thermal average, we use O(10 7 ) MC samplings. 
Hereafter, we will set J3/J1 = —0.1 in the model (|2j) in 
order to see effects of "small" additional interaction J3 
on defects of the J1-J2 model (pQ). 

Let us first see temperature ksT and anisotropy J2/J1 
dependences of spin and defect correlation functions 
along the x direction. These two correlators are defined 



G x (r x ) = (S(i >rj/ )5( rsB+ i >ri/ )) - (5(i >ri/ ))(5( ra , + i >ry ))(3) 

D x (r x ) = (^(1,^)^+1,^)) - (^(l,r y ))(^(r x +l,r y ))- (4) 

The quantity d r = (1 — S r S( rx ^ y +i^)/2 stands for the de- 
fect operator on the bond between (r Xl r y ) and (r x ,r y + 
1). It takes unity when a defect exists on the bond, while 
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FIG. 2: (a)(b) Spin correlation functions G x (r x ) of the mod- 
els (pQ) and J2}. (c)(d) Defect correlation functions D x (r x ) of 
the models (pQ) and ©. We have set L — 128 and plotted 
only the results of even r x . The amplitudes of odd r x are 
very small. 



it becomes zero otherwise. Figure [2] shows the T depen- 
dence of both spin and defect correlation functions in 
the models (pQ) and at J\/J 2 = 0.7. The panels (a) 
and (c) clearly indicates that spin and defect correlation 
functions decay exponentially irrespective of ksT in the 
model (pQ) without J3 term. This is consistent with the 
fact that the J1-J2 model has no phase transition. On the 
other hand, panels (b) and (d) show that the decay fash- 
ion of spin and defect correlators changes an algebraic 
type from an exponential form when T is sufficiently de- 
creased in the J1-J2-J3 model. This power-law behavior 
strongly suggests the emergence of a quasi long-range or- 
der of defects and spins, namely, a Berezinskii-Kosterlitz- 
Thouless (BKT) phase 0-13. 

To judge whether or not a BKT phase exists and to 
correctly determine the BKT transition point from the 
MC results, we should more carefully treat them rather 
than those of usual second-order transitions. Here, we 
utilize the finite-size scaling based on ratios of two spin 
correlators, that has been developed in Ref.[]jt Around a 
continuous transition or in a BKT phase of the models (pQ) 
and (p|) with size L (if they exist), a ratio of two spin 
correlation functions is expected to satisfy the following 
scaling relation: 



G x (r x ,T,L) 



f{L/£ x ), 



(5) 



where f(L/£ x ) is a function, £ x is the correlation length of 
G x , and the ratio of two lengths r x / r' x is fixed regardless 
of T and L. Here, we have explicitly denoted T and L 
dependences of G x . Since £ x diverges in a BKT phase, 
the left-hand side of Eq. (J5j) becomes independent of the 
size L there. In Fig. [3fa) and (b), we plot 



R{L) 



G X (L/2,T,L) 
G X (L/4,T, L) 



(6) 
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FIG. 3: (a)(b) Ratio of two spin correlators R(L) as a function 
of T. (c)(d) R(L) as a function of L/^ x with £ x 
Parameters a and t (i.e., Tbkt) are finely tuned. 
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FIG. 4: Phase diagram of the model %j 1 -j 2 -j 3 in the space 
(fcfiT, J2). Incommensurate ordered phases (see the text) 
might appear between BKT and AF ordered phases. 



as a function of T in the J1-J2-J3 models with different 
sizes. These two panels uncover that R(L) is almost in- 
variant under changing the size L in a finite temperature 
range. This strongly suggests that a BKT phase emerges 
due to the J3 term. In the paramagnetic phase near the 
BKT transition point, £ x is expected to behave as 



£ x oc exp(a/v / t) 



(7) 



where a is a nonuniversal constant, t = (T — 
Tbkt) /Tbkt-, and Tbkt is the BKT transition tempera- 
ture. Therefore, if we plot R(L) as a function of Le~ a ^ 
finely tuning parameters a and Tbkt-, R(L) with differ- 
ent L and T all align on a single curve in the regime 
Such curves are depicted in Fig.[3^c) and (d). Us- 
ing this scaling method, we determine Tbkt in the wide 
range of < J2/J1 < 1 with fixing J3/J1 = —0.1. These 
points Tbkt are given by filled circles in Fig. 2J 
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FIG. 5: Critical exponents of spin and defect correlation func- 
tions (a) ?7s and (b) rfd at J2/J1 = 0.7, where ^bTbkt / Ji ~ 
0.63. In the calculation of 77 S;C z, we have used the data G x (r x ) 
and D x (r x ) in the range 10 < r x < 50 for the L = 128 sys- 
tems. 



From spin and defect correlation functions in Fig. [2j 
we can also evaluate their critical exponents r] s and rjd, 
that are defined as 



G x (r x ) 



D x (r x ) 



(8) 



where spatially oscillating factors of G x and D x are ne- 
glected (We again note that all the plotted positions r x 
are even in Fig. [2]). The evaluated exponents are plotted 
in Fig. [5l Although there exist large errors in the data 
of Fig.[5l the result indicates that r] s (rjd) takes the value 
~ 1/4 (unity) at T = Tbkt-, and it monotonically de- 
creases with lowering T in the BKT phase. We remark 
that for the well-studied XY model on a square lattice, 
the critical exponent of the spin correlation also takes 1/4 
at the BKT transition point 15, 1(3, 18|. Namely, Ising 
spins S r of the present J1-J2-J3 model exhibit the same 
asymptotic behavior as XY spins of the XY model in the 
BKT phase. From this property of r] s ,d, we can also de- 
termine Tbkt as the point where r] s = 1/4 (rjd = 1). 
Temperatures Tbkt evaluated from 77^ are plotted by 
open circles in Fig. [H From these results, we conclude 
that a small J3 term causes a BKT phase, and it can 
be understood as the quasi long-range ordered phase of 
defects and spins. 

In addition to the BKT phase, we also find an AF or- 
dered phase in the J \- J 2- J 3 model. In fact, at an extreme 
case with J 2 = 0, the system is reduced to two decou- 
pled spatially anisotropic Ising models on square lattice 
that are exactly solvable and have a second-order phase 
transition to a Neel phase [20]. At J 2 = 0, the transition 
temperature T c is exactly calculated as 



sinh ( 

\kBT c 



S inh m 

KkbTc 



1, 



(9) 



which leads to k B T c w 0.906Ji at J 3 /Ji = -0.1. 
The ordering pattern is illustrated in Fig. [6] (cl) and 
(c2). Hence, their order parameters are respectively 
defined as mi = L~ 2 ^2 rx Ty (— 1) Ty S r and 777,2 = 
T~ 2 J2 r r (— l) ry+rx+1 5V" The corresponding suscepti- 
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FIG. 6: (a)(b) Susceptibility Xmax for the order parameter 
of the AF phase. (cl)(c2) AF ordering patterns. (dl)(d2) 
Examples of possible incommensurate ordering patterns. 
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FIG. 7: (a)(b) Spin structure factors S y (q y ). (c)(d) Temper- 
ature dependence of a peak position q P (< ty) of S y (q y ). 



bilities are given by 



Xl,2 



(K,2> - ( m ^f 



(10) 



In Fig. [6^a) and (b), we draw the MC evaluated sus- 
ceptibility Xmax = Max(xi,X2) as a function of T. It 
exhibits a peak at a certain temperature, and the peak 
becomes sharper with increasing the system size L. We 
have confirmed that below the temperature showing the 
peak, an AF order pattern (cl) or (c2) in Fig. [6] appears 
in the MC snap shots. Since the universality class of 
this transition is never known, we determine the tran- 
sition temperature T c through a naive finite-size scaling 
for the peak of Xmax- The evaluated T c are depicted in 
Fig. HI and one sees from it that T c is lower than Tbkt 
and it monotonically decreases with increasing J2/J1. 
When J2/J1 exceeds ~ 0.4, the evaluated T c becomes 
extremely small and the error of the finite-size scaling be- 
comes much larger. Therefore, we cannot determine T c 
in the range J2/J1 > 0.4. We also see that the peak posi- 
tions of Xmax seem to irregularly change as a function of 
L. This strange behavior of Xmax and large errors of the 
finite-size scaling for T c would be related with following 
things: (i) the MC simulation becomes less reliable in a 
low-temperature regime due to effects of frustration and 
defects, and (ii) the BKT phase above T c might largely 
affect the AF ordered phase when we consider finite-size 
systems. We will again discuss the nature of the AF tran- 
sition at T c later. We note that it is difficult to determine 
whether or not the BKT transition curve starts from the 
J2/J1 = axis and it coincides with T c at the axis in 
Fig.ffl 

Now, we turn to correlations along the y direction, es- 
pecially, focusing on the BKT phase. Similarly to Eq. (j4]), 
we define the spin correlation function along the y direc- 
tion G y (r y ) = (S(r x ,l)S(r Xj r y +l)) ~ ,1) ) ,r v +l) ) • 

Figure 0(a) and (b) show the spin structure factor for 



the y direction, that is defined as 
1 L_1 

Sy(q y ) = J COS fe r 2/) G 2/( r 2/)' ( U ) 

r y =0 

It has a double peak structure around wave number 
q y = 7T. The temperature dependence of one peak po- 
sition q y = q p < 7r are depicted in Fig. 0(c) and (d). 
They show that q p monotonically approaches to 7r with 
lowering T in the whole region < J2/J1 < 1- This 
is naturally understood because q p = tt in the AF or- 
dered phase (T < T c ), and q p is expected to deviate in 
proportion to the defect density (d r ) that monotonically 
grows with increasing T. From Fig. we see that G y is 
incommensurate in the paramagnetic and BKT phases. 
We have confirmed that the "commensurate" spin cor- 
relator G y (r y )/ cos(q p r y ) exhibits a power-law behavior 
~ ry Vs in the BKT phase. This is an additional evidence 
for the existence of the BKT phase. The incommensu- 
rate nature of G y also suggests the possibility of single 
or multiple intermediate phases with an incommensurate 
order between the BKT and AF ordered phases. Typical 
order patterns of possible intermediate phases are given 
in Fig. [6](dl) and (d2), and they are nothing but defect 
crystals. It is however hard to judge whether or not such 
phases exist since the MC simulation does not work well 
in ksT <C J\ and anyone never knows the universality of 
the phase transition to such intermediate phases. 

Viewpoints of statistical physics and field theory.— Fi- 
nally, we discuss the phases and transitions of the J±- 
J2-J3 model from the viewpoints of field theory and sta- 
tistical physics. We have uncovered that the J1-J2-J3 
model with J3/J1 = —0.1 has two successive transitions 
(Tbkt and T c ) with lowering T at least in the range 
0.1 < J2/J1 < 0.4. In addition to the present model, a 
few two-dimensional (2D) classical systems with a BKT 
phase in an intermediate-temperature range have been 
known. TV-state clock models with N > 5 are typical ex- 
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ample [27H29j| . Usually, in such systems, both the tran- 
sition between paramagnetic and BKT phases and that 
between BKT and ordered phases belong to the BKT 
universality class. Order parameters usually grows very 
slowly in the BKT class. However, the sharp suscep- 
tibility of Fig. [6] implies that the transition at T c is a 
first- or second-order type. To resolve this discrepancy, 
we can expect some scenarios. The first scenario is that 
the behavior of Xmax reduces to a BKT type in the ther- 
modynamic limit. The second is (as we mentioned) the 
possibility of additional incommensurate phases between 
BKT and AF phases. If the high-temperature transition 
is a BKT type, the low-temperature one may be first- 
or second-order. The third is that the transition between 
BKT and AF phases is a first-order type. In general, first- 
order transitions do not need to obey a usual picture of 
renormalization group. The fourth is that the AF order 
parameter and the defects are effectively decoupled 
in the long-distance physics. To determine the true sce- 
nario is an interesting future problem. To find the vortex 
picture of our BKT phase is another important issue. In 
fact, for the BKT (floating) phase of the 2D anisotropic 
next-nearest-neighbor Ising (ANNNI) model, the vortex 
has been defined For the BKT phase of the tri- 

angular Ising model with isotropic nearest-neighbor and 
next-nearest-neighbor couplings 24], [25| , a vortex picture 
has also been proposed [26|. 

BKT phases should be generally described by a confor- 
mal field theory (CFT) with central charge c = 1 30|, 3l| . 



The effective action for our BKT phase would be given 

by 



A = J dxdy ^(d^cp) 2 + gi cos(V4nK(f>) 



(12) 



where (x, y) is the continuous coordinate stemming from 
the discrete site index (r x ,r y ), <\> is a scalar field, and K 
is the so-called Tomonoga-Luttinger (TL) liquid param- 
eter [31]. In the BKT phase the g\ term with scaling 
dimension K and other perturbations are all irrelevant, 
and the effective action is given by the Gaussian part 
(9 M 0) 2 , while the g\ term becomes relevant [K < 2) in 
T > Tbkt and any quasi long-range correlation disap- 
pears. As we mentioned, Fig. [5] indicates that Ising spins 
play the same role as spins of 2D XY model in the BKT 
phase. The figure also suggests the relation rjd ~ 4t? s . 
From these two results, we can expect the following cor- 
respondence between operators of the model (j2j) and the 
c = 1 CFT: 



St* 



COS 



(y/ir/KO), d r ~ cos(2y/^/K0) (13) 



where 
tors. 



we have neglected 
The field is the 




spatially oscillating fac- 
dual |3fl, El of (/>, and 

r -n 2 /{2K) in the BKT 



Conclusion.— In conclusion, we have studied the spa- 
tially anisotropic triangular Ising model (j2j) with an ad- 
ditional Js term. We have shown that a small J3 < 
leads to a BKT phase in a wide parameter range and its 
effective field theory has been proposed. The BKT phase 
can be regarded as a quasi long-range ordered phase of 
defects and spins, and an incommensurate spin correla- 
tion along the y direction appears there. To develop the 
vortex picture and to reveal the low-temperature phases 
are important remaining issues. 
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phase. Since the low-temperature nature around T = T C 
has not been revealed enough, it is hard to develop the 
effective theory describing the physics around T c . 
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